#include "fermiqcd.h"

void main(int argc, char** argv) {
  string info=argv[1];
  mpi.open_wormholes(argc, argv);
  define_base_matrices("FERMILAB");

  int box[4];
  coefficients coeff;

  string gauge_filename=          prompt(info, "INPUT_GAUGE_FILENAME", "COLD");
  string fermi_filename=          prompt(info, "OUTPUT_FERMI_PREFIX", "fermi");

  box[0]=               (int) val(prompt(info, "T", "16"));
  box[1]=               (int) val(prompt(info, "X", "4"));
  box[2]=               (int) val(prompt(info, "Y", "4"));
  box[3]=               (int) val(prompt(info, "Z", "4"));

  string action_type=             prompt(info, "ACTION_TYPE", "CLOVER");
  string inverter=                prompt(info, "INVERTER", "MINRES");
  mdp_real absolute_precision=val(prompt(info, "ABSOLUTE_PRECISION", "0"));
  mdp_real relative_precision=val(prompt(info, "RELATIVE_PRECISION", "1e-6"));

  coeff["kappa_s"]=           val(prompt(info, "KAPPA_S", "0.13"));
  coeff["kappa_t"]=           val(prompt(info, "KAPPA_T", "0.13"));
  coeff["r_s"]=               val(prompt(info, "R_S", "1.0"));
  coeff["r_t"]=               val(prompt(info, "R_T", "1.0"));
  coeff["c_{sw}"]=            val(prompt(info, "C_SW", "0.0"));
  coeff["c_E"]=               val(prompt(info, "C_E", "1.0"));
  coeff["c_B"]=               val(prompt(info, "C_B", "1.0"));

  mpi << coeff;

  if(action_type=="CLOVER") 
    default_fermi_action=FermiCloverActionFast::mul_Q;
#if defined(SSE2)
  else if(action_type=="CLOVERSSE") 
    default_fermi_action=FermiCloverActionSSE2::mul_Q;
#endif
  else
    error("unknown action");

  if(inverter=="MINRES") 
    default_fermi_inverter=MinimumResidueInverter<fermi_field,gauge_field>;
  else if(inverter=="BICGSTAB") 
    default_fermi_inverter=BiConjugateGradientStabilizedInverter<fermi_field,gauge_field>;
  else
    error("unknown inverter");
		 
  int i,a,t;
  mpi << "Local Random Generator disable for space\n";
  mdp_lattice lattice(4,box,default_partitioning0,torus_topology,0,2,false);
  gauge_field U(lattice,3);
  fermi_field psi(lattice,3);
  fermi_field chi(lattice,3);
  fermi_field phi(lattice,3);
  site x(lattice);
  mdp_array<double,1> vector(box[0]);
  mdp_array<double,1> pseudoscalar(box[0]);

  if(gauge_filename=="HOT") set_hot(U);
  else if(gauge_filename=="COLD") set_cold(U);
  else U.load_as_float(gauge_filename);

  cout << "done\n";

  if(coeff.has_key("c_{sw}")) compute_em_field(U);

  cout << "done\n";

  for(t=0; t<box[0]; t++) vector(t)=pseudoscalar(t)=0;

  for(a=0; a<4; a++)
    for(i=0; i<3; i++) {
      psi=0;
      if(on_which_process(lattice,0,0,0,0)==ME) {
	x.set(0); 
	psi(x,a,i)=1;
      }
      psi.update();
      
      mul_invQ(chi,psi,U,coeff, absolute_precision,relative_precision);
      forallsites(x) 
	pseudoscalar(x(0))+=real(trace(chi(x)*hermitian(chi(x))));
      
      forallsites(x) 
	psi(x)=Gamma[1]*Gamma5*psi(x);
      mul_invQ(phi,psi,U,coeff,absolute_precision,relative_precision);
      forallsites(x)
	vector(x(0))+=real(trace(chi(x)*hermitian(Gamma5*Gamma[1]*phi(x))));
    }
  
  
  mpi.add(pseudoscalar.address(), box[0]);
  mpi.add(vector.address(), box[0]);

  FILE* fp;
  if(ME==0) {
    fp=fopen(fermi_filename.c_str(), "w");
    for(t=0; t<box[0]; t++)
      fprintf(fp, "%i,\t%e,\t%e\n", t, pseudoscalar(t), vector(t));
    fclose(fp);
  }
  mpi.close_wormholes();
}
